{ "cells": [ { "cell_type": "markdown", "id": "15f7c58f-7871-4db6-8861-63d422a7c5d5", "metadata": {}, "source": [ "The Community Aerosol and Radiation Model for Atmospheres, [CARMA](https://github.com/ESCOMP/CARMA), can be configured and run through musica. This tutorial walks through creating two small carma test..\n", "\n", "First, we need to import our packages" ] }, { "cell_type": "code", "execution_count": 1, "id": "8518f967-602a-450e-bf89-6f6b9037cc78", "metadata": {}, "outputs": [], "source": [ "import musica\n", "import ussa1976\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import xarray as xr\n", "from musica.constants import GAS_CONSTANT\n", "import musica.mechanism_configuration as mc\n", "import pandas as pd" ] }, { "cell_type": "markdown", "id": "310d4404-a3a0-45fb-949a-535f40683dc0", "metadata": {}, "source": [ "# Aluminum\n", "\n", "The first test is a coagulation test for aluminum in 5 bins. Please note that all input data going into the musica carma package should be in SI units, so meters instead of centimeters, pascal instead of hectapascals, etc. Anytime CARMA would be expecting a different unit, the musica package will automatically make that conversion when configuring carma.\n", "\n", "First up, we create the groups and elements, and create our carma paramters class. This is then used to make an instance of CARMA." ] }, { "cell_type": "code", "execution_count": 2, "id": "aa06c252-1b3b-4945-abec-f0ff97c0a45e", "metadata": {}, "outputs": [], "source": [ "group = musica.carma.CARMAGroupConfig(\n", " name=\"aluminum\",\n", " shortname=\"PRALUM\",\n", " rmrat=2.0,\n", " rmin=21.5e-6,\n", " rmon=21.5e-6,\n", " ishape=musica.carma.ParticleShape.SPHERE,\n", " eshape=1.0,\n", " mie_calculation_algorithm=musica.carma.MieCalculationAlgorithm.TOON_1981,\n", " is_ice=False,\n", " is_fractal=True,\n", " do_wetdep=False,\n", " do_drydep=True,\n", " do_vtran=True,\n", " solfac=0.0,\n", " scavcoef=0.0,\n", " df=[1.6] * 5, # 5 bins with fractal dimension 1.6\n", " falpha=1.0\n", ")\n", "\n", "# Create aluminum element\n", "element = musica.carma.CARMAElementConfig(\n", " igroup=1,\n", " isolute=0,\n", " name=\"Aluminum\",\n", " shortname=\"ALUM\",\n", " itype=musica.carma.ParticleType.INVOLATILE,\n", " icomposition=musica.carma.ParticleComposition.ALUMINUM,\n", " rho=0.00395, # kg/m3\n", " arat=[1.0] * 5, # 5 bins with area ratio 1.0\n", " kappa=0.0,\n", ")\n", "\n", "# Create coagulation\n", "coagulation = musica.carma.CARMACoagulationConfig(\n", " igroup1=1,\n", " igroup2=1,\n", " igroup3=1,\n", " algorithm=musica.carma.ParticleCollectionAlgorithm.FUCHS)\n", "\n", "params = musica.carma.CARMAParameters(\n", " nbin=5,\n", " nz=1,\n", " dtime=1800.0,\n", " groups=[group],\n", " elements=[element],\n", " coagulations=[coagulation]\n", ")\n", "\n", "FIVE_DAYS_IN_SECONDS = 432000\n", "params.nstep = FIVE_DAYS_IN_SECONDS // params.dtime\n", "params.initialization.do_vtran = False\n", "\n", "carma = musica.CARMA(params)" ] }, { "cell_type": "markdown", "id": "fd1fc77d-2bbd-4257-a480-6e7543f7bb69", "metadata": {}, "source": [ "Next, we make some envionrmental data which we will use to initialize our state." ] }, { "cell_type": "code", "execution_count": 3, "id": "1ca883d3-c507-43d0-8aba-e5ff8ea84312", "metadata": {}, "outputs": [], "source": [ "n_levels = params.nz\n", "deltaz = 1000.0\n", "zmin = 16500.0\n", "\n", "vertical_center = zmin + (np.arange(n_levels) + 0.5) * deltaz\n", "vertical_levels = zmin + np.arange(n_levels + 1) * deltaz\n", "\n", "centered_variables = ussa1976.compute(z=vertical_center, variables=[\"t\", \"p\", \"rho\"])\n", "edge_variables = ussa1976.compute(z=vertical_levels, variables=[\"p\"])\n", "\n", "temperature = centered_variables.t.values\n", "pressure = centered_variables.p.values\n", "pressure_levels = edge_variables.p.values\n", "density = centered_variables.rho.values" ] }, { "cell_type": "markdown", "id": "e9bb3e37-78cb-4831-8d4e-b4d677c22ae9", "metadata": {}, "source": [ "Finally, we will run the carma box model. While the model runs, we will also collect the data. At the end we will combine it all into one xarray dataset." ] }, { "cell_type": "code", "execution_count": 4, "id": "8d77d442-68ed-4640-8992-8f96ab823975", "metadata": {}, "outputs": [], "source": [ "mmr_initial = 5e9 / (deltaz * 2.57474699e14) / density[0]\n", "\n", "state = carma.create_state(\n", " time_step=params.dtime,\n", " temperature=temperature,\n", " pressure=pressure,\n", " pressure_levels=pressure_levels,\n", " vertical_center=vertical_center,\n", " vertical_levels=vertical_levels,\n", " longitude=0.0,\n", " latitude=-105.0,\n", " coordinates=musica.carma.CarmaCoordinates.CARTESIAN,\n", ")\n", "\n", "for i in range(params.nbin):\n", " for j in range(len(params.elements)):\n", " state.set_bin(i + 1, j + 1, mmr_initial)\n", "\n", "bin_data = state.get_bins()\n", "bin_data = bin_data.expand_dims({\"time\": [0]})\n", "env = state.get_environmental_values()\n", "env = env.expand_dims({\"time\": [0]})\n", "time_array = [0.0] # Start with time 0\n", "\n", "# Run the simulation for the specified number of steps\n", "for step in range(1, int(params.nstep)):\n", " state.step()\n", " bin_data = xr.concat([bin_data, state.get_bins().expand_dims({\"time\": [step * params.dtime]})], dim=\"time\")\n", " env = xr.concat([env, state.get_environmental_values().expand_dims(\n", " {\"time\": [step * params.dtime]})], dim=\"time\")\n", " time_array.append(step * params.dtime)\n", "\n", "ds = xr.merge([bin_data, env])" ] }, { "cell_type": "code", "execution_count": 5, "id": "3790977e-ac2c-4289-95e1-cdc5c13c834e", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.Dataset> Size: 154kB\n", "Dimensions: (time: 240, bin: 5, element: 1,\n", " vertical_center: 1, vertical_level: 2)\n", "Coordinates:\n", " * time (time) float64 2kB 0.0 1.8e+03 ... 4.302e+05\n", " * bin (bin) int64 40B 1 2 3 4 5\n", " * element (element) int64 8B 1\n", " * vertical_center (vertical_center) float64 8B 1.7e+04\n", " * vertical_level (vertical_level) float64 16B 1.65e+04 1.75e+04\n", "Data variables: (12/18)\n", " mass_mixing_ratio (time, bin, element, vertical_center) float64 10kB ...\n", " number_mixing_ratio (time, bin, element, vertical_center) float64 10kB ...\n", " number_density (time, bin, element, vertical_center) float64 10kB ...\n", " nucleation_rate (time, bin, element, vertical_center) float64 10kB ...\n", " wet_particle_radius (time, bin, element, vertical_center) float64 10kB ...\n", " wet_particle_density (time, bin, element, vertical_center) float64 10kB ...\n", " ... ...\n", " sedimentation_flux (time, bin, element) float64 10kB 0.0 ... 0.0\n", " deposition_velocity (time, bin, element) float64 10kB -9.99 ... -...\n", " temperature (time, vertical_center) float64 2kB 216.6 ......\n", " pressure (time, vertical_center) float64 2kB 8.85e+03 ...\n", " air_density (time, vertical_center) float64 2kB 0.1417 .....\n", " latent_heat (time, vertical_center) object 2kB None ... None
<xarray.Dataset> Size: 177kB\n", "Dimensions: (time: 240, bin: 5, element: 1,\n", " vertical_center: 1, vertical_level: 2,\n", " gas: 2)\n", "Coordinates:\n", " * time (time) float64 2kB 0.0 1.8e+03 ... 4.302e+05\n", " * bin (bin) int64 40B 1 2 3 4 5\n", " * element (element) int64 8B 1\n", " * vertical_center (vertical_center) float64 8B 1.7e+04\n", " * vertical_level (vertical_level) float64 16B 1.65e+04 1.7...\n", " * gas (gas) <U5 40B 'H2O' 'H2SO4'\n", "Data variables: (12/24)\n", " mass_mixing_ratio (time, bin, element, vertical_center) float64 10kB ...\n", " number_mixing_ratio (time, bin, element, vertical_center) float64 10kB ...\n", " number_density (time, bin, element, vertical_center) float64 10kB ...\n", " nucleation_rate (time, bin, element, vertical_center) float64 10kB ...\n", " wet_particle_radius (time, bin, element, vertical_center) float64 10kB ...\n", " wet_particle_density (time, bin, element, vertical_center) float64 10kB ...\n", " ... ...\n", " gas_mass_mixing_ratio (time, gas, vertical_center) float64 4kB ...\n", " gas_saturation_wrt_ice (time, gas, vertical_center) float64 4kB ...\n", " gas_saturation_wrt_liquid (time, gas, vertical_center) float64 4kB ...\n", " gas_vapor_pressure_wrt_ice (time, gas, vertical_center) float64 4kB ...\n", " gas_vapor_pressure_wrt_liquid (time, gas, vertical_center) float64 4kB ...\n", " weight_pct_aerosol_composition (time, gas, vertical_center) float64 4kB ...